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SUMMARY 


This  report  presents  and  studies  a computational  method  called  "Glimm's 
method"  from  the  point  of  view  of  applications  to  blast  wave  problems.  The 
report  includes  the  results  of  preliminary  numerical  studies.  This  work  was 
performed  by  the  author  who  is  a member  of  the  Applied  Mathematics  Branch  of 
NAVSWC/WOL.  Conversations  with  Professor  Alexandre  Chorin,  Professor  Gary  Sod, 
Dr.  Hy  Sternberg,  and  Dr.  Greg  Shubin  are  gratefully  acknowledged. 

PAUL  R.  WESSEL 
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INTRODUCTION 


Glimm's  method  is  a numerical  technique  for  solving  hyperbolic  systems  of 
conservation  laws  in  one  space  dimension.  It  has  been  extended  to  problems  with 
more  than  one  space  dimension  and  to  systems  in  weak  conservation  form.  The  method 
does  not  involve  finite  differences;  its  unique  features  are  a random  choice  at 
each  time  step  (which  implies  that  the  solution  is  a random  variable)  and  the 
explicit  solution  of  a Riemann  (shock-tube)  problem  for  each  pair  of  adjacent  mesh 
points.  The  main  advantage  of  the  method  is  its  ability  to  resolve  discontinuities 
(both  shock  waves  and  contact  discontinuities)  exactly  without  smoothing.  Further- 
more, this  takes  place  automatically — no  special  shock  tracking  procedures  are 
required.  In  particular,  the  method  recognizes  the  formation  of  a discontinuity 
and  correctly  resolves  the  interaction  of  discontinuities.  The  price  that  the  user 
pays  for  this  property  is  the  imposition  of  a statistical  error  in  the  computed 
solution. 

A precursor  of  Glimm's  method  is  the  finite  difference  scheme  of  Godunov. 
Riemann  problems  are  solved  as  in  the  Glimm  algorithm;  however,  there  is  no  random 
choice  element  and  the  solution  is  advanced  in  time  by  a deterministic  scheme 
[Reference  1].  The  complete  algorithm  was  introduced  by  Glimm  [Reference  2]. 

Glimm  used  the  algorithm  as  a theoretical  tool  to  be  used  in  obtaining  existence 
(for  all  time)  proofs  for  (weak)  solutions  of  hyperbolic  systems  of  conservation 
laws.  As  he  defined  the  algorithm,  the  convergence  rate  is  too  slow  for  numerical 
computations;  nevertheless,  the  paper  represents  a major  theoretical  advance.  This 
state  of  affairs  did  not  change  substantially  until  Chorin's  work  was  published  in 
1976  [Reference  3].  In  this  paper,  Glimm's  algorithm  was  modified  in  such  a way 
that  it  became  acceptable  from  the  point  of  view  of  computational  efficiency,  and 
yet  its  theoretical  properties  remained  unchanged. 

This  report  will  not  consider  the  theoretical  aspects  of  the  theory  of 
hyperbolic  conservation  laws,  the  existence  and  uniqueness  of  solutions  to  Riemann 
problems,  or  the  convergence  properties  of  Glimm's  method  applied  to  such  conserva- 
tion laws.  There  is  an  extensive  literature  on  these  topics  and  we  refer  the 
reader  to  two  papers  of  Lax  [References  4 and  5];  the  later  paper  contains  a 
substantial  bibliography. 

1.  Godunov,  S.  K. , "Finite  Difference  Methods  for  Numerical  Computation  of 
Discontinuous  Solutions  of  the  Equations  of  Fluid  Dynamics,"  Mat.  Sbornik  47, 
1959,  p.  271. 

2.  Glimm,  J.,  "Solutions  in  the  Large  for  Nonlinear  Hyperbolic  Systems  of  Equations," 
C.P.A.M.  18.  1965,  p.  697. 

3.  Chorin,  A.  J.,  "Random  Choice  Solution  of  Hyperbolic  Systems,"  J.C.P.  22, 

Dec.  1976,  p.  517-533. 

4.  Lax,  P.  D.,  "Hyperbolic  Systems  of  Conservation  Laws  II,"  C.P.A.M.  10,  1957, 
p.  537-566. 

5.  Lax,  P.  D.,  Hyperbolic  Systems  of  Conservation  Laws  and  the  Mathematical  Theory 
of  Shock  Waves.  Society  for  Industrial  and  Applied  Mathematics,  1973. 
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Our  main  objective  is  to  indicate  how  Glimm's  method  can  be  applied  to  blast 
wave  problems.  The  advantages  given  above  for  Glimm's  method  indicate  why  such 
an  application  would  be  desirable.  In  a blast  wave,  it  is  the  discontinuous 
portions  of  the  flow  field  which  are  of  the  greatest  importance;  any  method  specif- 
ically designed  for  these  regions  should  be  quite  useful.  Furthermore,  the 
statistical  error  which  is  inherent  in  the  method  will  not  be  important  in  most 
applications.  These  points  will  again  be  taken  up  in  the  last  section  of  this 
report. 

Blast  wave  problems  can  be  formulated  as  systems  of  hyperbolic  conservation 
laws,  either  strong  or  weak  depending  on  the  geometry  of  the  problem.  We  refer, 
of  course,  to  the  conservation  laws  of  gas  dynamics-conservation  of  mass,  momentum, 
and  energy.  Thus,  the  equations  of  gas  dynamics  will  be  the  only  application  of 
Glimm's  method  to  be  considered  here  in  any  detail.  In  the  next  section,  the 
basic  Glimm  algorithm  for  the  one-dimensional  planar  case  is  described.  Our 
treatment  is  based  on  two  papers  of  Chorin  [References  3 and  6].  However,  we 
include  a discussion  of  some  very  simple  examples  to  motivate  the  technique.  The 
next  section  takes  up  the  two-dimensional  planar  case  and  follows  Chorin  [Reference 
3].  The  equations  of  gas  dynamics  are  in  strong  conservation  form  for  planar 
symmetry.  This  is  no  longer  the  case  for  spherical,  cylindrical,  and  axisymmetric 
symmetry  and  the  extension  of  Glimm's  method  to  the  weak  conservation  form  case  is 
taken  up  in  the  following  section.  Throughout  these  three  sections,  the  results 
of  numerical  tests  will  be  presented;  this  will  include  some  new  results.  In  the 
final  section,  some  blast  wave  problems  will  be  discussed  from  the  point  of  view  of 
Glimm's  method.  One  interesting  point  in  such  applications  is  that  the  equations 
of  state  which  are  of  interest  can  be  quite  complicated.  However,  Glimm's  method 
has  thus  far  only  been  applied  to  a y-law  gas.  We  hope  that  this  gap  can  be 
considered  partially  filled  by  the  appendix  to  this  report  where  we  show  how 
Riemann  problems  can  be  solved  for  an  arbitrary  equation  of  state. 

GLIMM'S  METHOD 

In  this  section,  the  Glimm  algorithm  for  solving  initial  value  problems  for 
strong  conservation  laws  in  one  space  dimension  is  presented.  The  problem  to  be 
solved  is 


u + f(u)  = 0 
— t — — x 

u(x,0)  given 


(1) 


where  u * u(x,t)  is  a function  defined  on  the  half-plane  t 1 0.  Physically,  the 
vector  ii  represents  the  conserved  quantities  and  the  vector  jf  represents  the 
fluxes.  The  system  (1)  is  the  differential  version  of  the  integral  conservation 
laws.  The  equations  of  one-dimensional  gas  dynamics  with  planar  symmetry  are 


f 


6.  Chorin,  A.  J.,  "Random  Choice  Methods  with  Applications  to  Reacting  Gas  Flow," 
J.C.P.  25,  Nov.  1977,  p.  253-272. 
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Pt  + (pu)^  = 0 


(pu)t  + (pu  + p>x  * 0 


et  + ((e  + P)u)x  = 0 


e = pc  + %pu  , e = F (p,p) 


where  p = density,  u = velocity,  e = energy  per  unit  volume,  e * internal  energy 
per  unit  mass,  and  p = pressure.  The  equation  e = F(p,p)  is  the  equation  of  state 
and  we  assume  that  this  relation  can  be  inverted  to  find  the  pressure  p as  a 
function  of  e and  p.  If  we  set  u = (p,pu,e)T,  it  is  clear  that  the  system  (2)  can 
be  put  in  the  form  (1)  after  using  the  equation  of  state  to  eliminate  the  pressure 
terms.  A computation  shows  that  the  system  (2)  is  hyperbolic  for  any  reasonable 
equation  of  state. 


The  method  of  computation  for  solving  the  system  (1)  is  as  follows.  Let  k be 
an  increment  of  time  and  h a spatial  increment.  The  solution  is  to  be  obtained  at 
the  times  t = nk  at  the  points  x = ih,  i = 0,±1,±2,***  and  at  the  times  t = (n+^k 
at  the  points  x = (i+}5)h,  i = 0,±1,±2,**»  is  here  n ranges  over  the  nonnegative 

integers.  Let  u,  ^ u(ih,nk),  uV. , 'v  u((i+^)h,  (n+%)k).  To  define  the  algorithm, 

— 1 i , i 

n+^  ri  n 

it  is  necessary  and  sufficient  to  describe  how  to  find  u^+^  given  u^,  u^+^. 
Consider  the  initial  value  problem  for  the  system  (1)  with  initial  data 


^ u((i+^)h,  (n+^k).  To  define  the  algorithm, 
j .f  i i * r -i i . . n+  n n 


for  x ^ 0 


u(x,0) 


for  x < 0. 


This  is  called  a Riemann  problem.  Let  v(x,t)  denote  the  solution  of  this  problem 
and  let  9^  be  a value  of  a random  variable  9 equidistributed  in  [-h,h] • Let 
H = v(6ih,k/2)  = the  value  of  the  solution  of  the  Riemann  problem  at  the  point 
(0ih,k/2).  Set 


n+h 


This  procedure  is  repeated  for  each  i and  a similar  construction  takes  place  from 
t * (n+Js)k  to  t * (n+l)k. 


This  construction  is  illustrated  in  Figures  1,  2,  and  3.  Figure  (la)  is  a 
schematic  diagram  of  a shock  tube  at  t * 0 with  constant  states  ju£,ur  to  the  left 
and  right  of  a membrane  located  at  x = 0.  Figure  (lb)  shows  the  same  shock  tube 
at  some  later  time  t * t*  after  the  membrane  has  been  burst.  In  this  illustration, 
a shock  is  moving  to  the  right  and  a rarefaction  is  moving  to  the  left.  Figure  (lc) 
shows  the  solution  in  the  (x.t)-plane;  note  that  the  solution  depends  only  on  the 
ratio  x/t.  Of  course.  Figure  1 is  only  one  example  of  a solution  to  a Riemann 
problem;  which  waves  form,  their  velocities,  and  the  values  of  the  state  variables 
in  the  *-state  must  be  determined  by  some  means.  Figure  2 illustrates  the 
situation  in  the  (p,u)-plane.  There  is  a one-parameter  family  of  states 


(p,u)  which  can  be  connected  to  each  of  (p  ,u.),  (p  ,u  ) 

k k r r 
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(1c)  THE  SOLUTION  IN  THE  lx,  t)  PLANE 


FIGURE  1.  SOLUTION  OF  THE  SHOCK  TUBE  PROBLEM  FOR  THE  CASE  OF  A 
FORWARD  SHOCK  AND  A BACKWARDS  RAREFACTION  . 
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by  either  a uniform  shock  or  a simple  wave  (centered  rarefaction  fan).  A general 
proof  of  this  fact  may  be  found  in  Reference  5.  These  one-parameter  families  are 
represented  by  the  curves  in  Figure  2.  The  intersection  of  the  two  curves  through 
(p^.uj)  and  (pr,ur)  must  be  (p*,u*);  once  the  latter  is  known,  the  complete 
Riemann  problem  solution  is  easy  to  find.  The  analytical  method  used  to  find 
(p*,u*)  - the  Godunov  iteration  - is  explained  in  detail  in  the  appendix  to  this 
report  after  which  the  complete  solution  is  constructed.  Most  of  the  ideas 
involved  are  explained  and  very  well  motivated  in  Reference  7.  The  theory  for 
general  hyperbolic  systems  of  conservation  laws  may  be  found  in  References  4 and  5. 


Assume  now  that  Riemann  problems  can  be  solved  for  the  system  (2).  The 
mechanics  of  Glimm's  method  are  illustrated  in  Figure  3.  At  time  t = nAt,  the 
Riemann  problem,  setting  the  values  of  the  state  variables  at  x = iAx  and 
x = (i+l)Ax  to  be  uj.Uj-,  is  constructed  with  the  membrane  at  x = (i+V)Ax.  The 
solution  to  this  Riemann  problem  is  found  (see  the  appendix)  and  is  then  sampled 
along  the  line  t = (n+^)At  where  the  sample  point  has  the  uniform  distribution  on 
the  segment  [ iAx, (i+1) Ax] . The  values  of  u at  the  sample  point  are  then  assigned 


n+^j 

-i+V 


This  procedure  is  then  repeated  for  each  pair  of  adjacent  mesh  points 


along  the  x-axis;  then  a new  random  number  is  selected  and  one  proceeds  to  the  next 
time  step.  For  a more  thorough  discussion  of  the  algorithm  and  its  properties, 
see  References  3 and  6. 


Glimm's  method  is  unconditionally  stable.  Nevertheless,  the  Courant  condition 
must  be  satisfied;  if  not,  waves  will  leave  the  sampling  interval  which  implies 
that  the  sampling  probabilities  will  be  incorrect.  In  effect,  this  changes  the 
problem  being  solved.  The  success  of  the  method  is  dependent  on  the  choice  of 
random  numbers  {0n).  Chorin  introduced  two  new  features  here — first,  only  one 
value  of  0 for  each  time  step  (instead  of  choosing  a new  random  number  for  each 
mesh  point  every  time  step  as  Glimm  did  in  his  paper)  is  used  and,  second,  special 
techniques  are  used  to  insure  that  the  sequence  approaches  rapidly  the  uniform 
distribution.  A recent  result  of  Liu  [Reference  8]  has  shown  that  it  is  not 
necessary  to  use  random  numbers  at  all;  all  that  is  necessary  is  that  the  sequence 
used  approach  equidistribution.  Thus,  in  the  near  future,  it  will  be  possible  to 
use  predetermined  sequences  which  are  chosen  to  optimize  the  approach  to  equi- 
distribution; this  development  will  improve  the  solutions  obtained  with  Glimm's 
method . 


Glimm's  method  is  only  first-order  accurate  but  has  infinite  resolution  in  a 
sense  to  be  described  by  a few  simple  examples.  In  Figure  (4a),  we  have  chosen 
two  constant  states  ti^  and  ur  in  such  a way  that  the  solution  to  the  associated 
Riemann  problem  is  a unit  speed  shock  moving  to  the  right  (with  no  wave  at  all 
moving  to  the  left  and  no  contact  discontinuity).  This  is  an  easy  case  to  analyze. 


7.  Courant,  R.  B.  , Friedrichs,  K.  0.,  Supersonic  Flow  and  Shock  Waves,  ;.’iley- 
Interscience,  1948. 

8.  Liu,  T.  P.,  "The  Deterministic  Version  of  the  Glimm  Scheme,"  Comm.  Math.  Phys. 
57,  1977,  p.  135-148. 
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(4a)  THE  CONFIGURATION  FOR  A REIMANN  PROBLEM  WHOSE  SOLUTION  IS  A SHOCK: 

THIS  IS  A LIMITING  CASE  IN  WHICH  THE  LEFT  WAVE  AND  CONTACT  DISCONTINUITY 
ARE  NOT  PRESENT. 


(4b)  LEFT,  LEFT  CHOICE 


* * * * * 

\\\i  w 

///!  // 

* * * ! * * « 

(4c)  LEFT,  RIGHT  CHOICE 


* * * 1 * ♦ * 

//■/// 

\\  KV\ 

* * * | * * * 

I4d)  RIGHT,  LEFT  CHOICE 


* * * * * * 


\ \ \i  \ \ 

* * * * * 

\ \ \ \ \ 

* * | * # * 

(4e)  RIGHT,  RIGHT  CHOICE 


FIGURE  4.  THE  SOLUTION  FOR  A UNIT  SPEED  UNIFORM  SHOCK  WAVE  USING  GLIMM'S  METHOD. 
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When  choosing  the  random  number  in  Gllmm's  method,  there  are  only  two  possibilities 
— either  u^  is  obtained  or  ur  is  obtained.  These  occur  precisely  when  the  random 
choice  lies  to  the  left  or  the  right  of  the  shock,  respectively.  For  two 
consecutive  half-time  steps,  there  are  four  possibilities  and  these  are  illustrated 
in  Figures  (4b)  - (4d).  For  example.  Figure  (4b)  illustrates  the  fact  that  if  two 
consecutive  random  numbers  are  chosen  to  the  left  of  the  shock,  then  the  shock 
moves  forward  one  mesh  point.  From  Figure  (4a)  it  is  clear  that  pr(uj^)  = probabil- 
ity that  the  random  number  lies  to  the  left  of  the  shock  « \ and 

At-  2Ax 

pr(V  = * ~ IS*  HenCe* 

Pr(ue)2  - Pr(ur)2  - 

This  says  that  the  net  probability  of  the  shock  advancing  one  mesh  point  per  full 
time  step  is  At /Ax.  The  velocity  inherent  in  an  advance  of  one  mesh  point  is  just 
Ax/At.  Therefore,  the 


mean  shock  speed  = (At /Ax) • (Ax/ At)  = 1. 

Hence,  on  the  average  the  shock  moves  at  the  correct  speed.  At  least  as  important 
as  this  result  is  the  fact  that  the  shock  remains  a shock  without  smearing  and  the 
values  of  the  flow  field  on  either  side  of  the  shock  remain  exact  (up  to  the 
tolerance  used  in  solving  the  Riemann  problems).  I.e.,  the  resolution  of  the 
method  is  infinite  for  this  example.  Even  though  the  location  of  the  shock  will 
generally  be  in  error,  one  will  always  see  a shock,  it  will  be  the  correct  shock, 
and  no  new  unwanted  waves  of  any  kind  will  be  introduced. 

This  argument  involving  the  shock  and  rarefaction  curves  in  the  (p,u)-plane 
can  be  generalized  to  show  that  any  Riemann  problem  will  be  solved  by  Glimm's 
method  with  infinite  resolution  despite  inherent  statistical  errors.  That  is, 
the  information  content  of  the  waves  is  never  lost. 

A somewhat  more  difficult  problem  is  the  interaction  of  the  waves  from  two 
adjacent  Riemann  problems.  The  setup  and  analytical  solution  of  a problem  of  this 
type  is  illustrated  in  Figure  5.  The  various  waves  remain  uniform  throughout  the 
solution.  We  have  also  solved  this  problem  using  Glimm's  method  (for  the  treatment 
of  the  boundaries,  which  is  straightforward  here,  see  Reference  3).  As  in  the 
previous  problem,  the  wave  speeds  are  subject  to  statistical  error  and  yet  the 
constant  states  are  computed  exactly  to  the  precision  of  the  Riemann  problem 
solver.  The  interaction  zone  just  below  state  E contained  only  a few  mesh 
points;  still  state  E was  computed  exactly.  This  illustrates,  in  this  more 
complicated  case,  the  infinite  resolution  of  the  algorithm  and  the  property  that 
it  does  not  lose  information. 

Of  course,  if  one  is  faced  with  arbitrary  initial  data,  the  solution  will 
rarely  consist  of  uniform  waves.  In  view  of  the  fact  that  the  basic  idea  of 
Glimm's  method  is  to  approximate  the  flow  field  by  a collection  of  simple  waves, 
the  question  arises  as  to  how  useful  the  method  will  be  for  strong  nonuniform 
waves.  Before  proceeding  to  a numerical  experiment,  we  comment  that  as  the  mesh 
spacing  becomes  smaller,  the  approximation  of  the  flow  field  by  piecewise  constant 
data,  hence  the  simple  wave  approximation  to  the  solution,  becomes  better.  Thus, 
we  expect  more  accurate  solutions  with  more  mesh  points — just  as  with  any  finite 
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FIGURE  5.  AN  EXAMPLE  OF  THE  INTERACTION  OF  UNIFORM  WAVES.  THE  CONS- 
TANT STATES  A AND  B ARE  THE  INITIAL  DATA;  THE  PICTURE  IS  SYM- 
METRIC ABOUT  THE  CENTER  OF  STATE  B THE  STATES  C AND  D ARISE 
FROM  THE  SOLUTION  OF  THE  RIEMANN  PROBLEM  GIVEN  BY  A AND  B. 
A'  IS  A CONSTANT  STATE  GIVEN  BY  THE  SHOCK  REFLECTION  EQUA- 
TIONS. STATE  E.  WHICH  ARISES  FROM  THE  INTERACTION  OF  THE 
RAREFACTION  FANS,  IS  ALSO  A CONSTANT  STATE.  FOR  AN  IDEAL  GAS 
WITH  7 = 3,  THE  SOLUTION  MAY  BE  OBTAINED  ANALYTICALLY. 
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difference  scheme.  Still,  this  was  not  the  case  In  the  above  simple  examples 
where  no  Information  was  lost  so  long  as  there  were  as  few  as  two  or  three  mesh 
points  for  each  distinct  wave  present  in  the  solution.  The  question  of 
convergence  and  convergence  rates  for  Glimm’s  method  in  such  problems  are 
unanswered;  neither  Glimm's  original  proof  nor  any  of  the  later  theoretical  work 
on  the  method  apply  to  the  case  of  gas  dynamics  with  general  initial  data.  Thus, 
we  pursue  a numerical  example. 

Fortunately,  exact  analytical  solutions  are  available  for  purposes  of  testing 
the  method.  We  have  pursued  an  example  of  an  exponential  shock  path  pushed  by  a 
piston  due  to  Sternberg.  The  wave  pattern  between  the  shock  and  the  piston  is 
about  as  nonuniform  as  possible.  The  exact  solution  may  be  found  in  Reference  9. 

A comparison  of  this  exact  solution  with  a computation  using  Glimm's  method  is 
illustrated  in  Figure  6.  Observe  that  both  the  piston  and  the  shock  have  moved 
correctly  and  that  the  wave  profile  is  correct.  We  hope  to  do  further  work  with 
similarity  solutions  along  these  lines,  especially  in  the  reactive  flow  case  (see 
Reference  6 for  an  extension  of  Glimm's  method  to  this  case;  also,  the  operator 
splitting  idea  to  be  presented  in  a later  section  can  be  applied).  Two  papers  of 
Sternberg  [References  9 and  10 1 will  be  the  basis  for  these  computations. 

This  completes  our  survey  of  the  basic  Glimm's  method  for  strong  hyperbolic 
conservation  laws  in  one  space  dimension.  As  a final  remark,  we  note  that  both 
the  theory  and  the  practice  of  the  method  are  easier  to  understand  for  single 
conservation  laws  than  for  systems.  Burger's  equation,  ut  + uux  ■ 0 is  a good 
example.  The  situation  considered  in  Reference  11  is  especially  interesting  in 
this  regard. 

GLIMM’S  METHOD  - EXTENSION  TO  TWO  INDEPENDENT  VARIABLES.  The  extension  of 
Glimm's  method  to  systems  of  strong  conservation  laws  in  more  than  one  space 
dimension  would  be  straightforward  if  multidimensional  Riemann  problems  could  be 
solved.  For  example,  the  two-dimensional  Riemann  problem  consists  of  the 
specification  of  four  constant  states,  one  for  each  quadrant  of  the  plane,  at 
t *■  0;  and  the  solution  to  this  problem  would  consist  of  using  this  configuration 
as  the  initial  state  and  finding  the  complete  solution  of  the  equations  in  a 
neighborhood  of  the  singularity  at  the  origin  for  all  t > 0.  Unfortunately,  very 
little  is  known  about  this  problem.  Also,  if  an  exact  solution  (or  iteration 
procedure)  is  ever  discovered,  it  is  bound  to  be  considerably  more  difficult  than 
the  one-dimensional  case. 

The  procedure  developed  by  Chorin  [Reference  3]  is  to  use  the  one-dimensional 
Glimm  algorithm  as  a building  block  in  a fractional  step  method.  In  place  of  the 
two  half-steps  of  the  one-dimensional  method,  one  takes  four  quarter  steps  for  the 
two-dimensional  method;  each  quarter  step  is  a sweep  in  either  the  x or  y direction 


9.  Sternberg,  H.  M. , "Similarity  Solutions  for  Reactive  Shock  Waves,"  Quart . 

J.  of  Mech.  and  Appl.  Math.  23.  Feb.  1970,  p.  77-99. 

10.  Sternberg,  H.  M. , "Constant  Velocity  Reactive  Shock  Waves  for  Testing 
Numerical  Methods,"  preprint,  1978. 

11.  Concus,  P.  and  Proskurowski,  W. , "Numerical  Solution  of  a Nonlinear  Hyperbolic 
Equation  by  the  Random  Choice  Method,"  Lawrence  Berkeley  Laboratory  Report 
LBL-6487  Rev.,  Dec.  1977. 
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This  idea  can  be  generalized  to  more  dimensions,  but  the  discussion  here  will  be 
restricted  to  the  two-dimensional  case.  The  equations  of  two-dimensional 
compressible  flow  with  planar  symmetry  are 

p + (pu)  + (pv)  - 0 
l x y 

(pu)  + (pu2  + p)  + (puv)  - 0 

t x y 

(pv)  + (puv)  + (pv2  + p)  - 0 (3) 

t x y 

e + ((e  + p)u)  + ((e  + p)v)  ■ 0 
t x y 

2 2 

e = PE  + ^p(u  + v ),  E - F(p,p). 

This  is  clearly  of  the  form 


“t  + I<u)x  + £(u)y  - 0 

T 

where  the  conserved  quantities  are  u = (p,pu,pv,e)  . As  before,  it  is  assumed  that 

the  equation  of  state  may  be  inverted  to  find  p as  a function  of  e and  p,  hence  as 

a function  of  the  conserved  quantities.  In  this  manner  the  flux  vectors  _f  and  £ 

may  be  found  from  (3).  The  equations  to  be  solved  in  an  x-sweep  are 

Pfc  + (pu)  “ 0 

(pu)  + (pu2  + p)  =0 

X (4) 

(pv)  ^ + (puv)x  = 0 

et  + ((e  + P>u)x  = 0. 

That  is,  all  y-derivatives  are  set  equal  to  zero.  The  third  of  these  equations 
becomes 


v + (uv)  = 0 (5) 

t x 

after  application  of  the  first  equation.  This  means  that  the  y-component  of 
velocity  v is  transported  as  a passive  scalar  in  an  x-sweep.  This  is  a particularly 
simple  situation  for  Glimm's  method;  the  solution  of  a Riemann  problem  for  the 
1st,  2nd,  and  4th  of  equation  (4)  is  independent  of  equation  (5).  Hence,  an 
x-sweep  reduces  to  the  system  (2)  of  one-dimensional  flow.  In  view  of  the  third 
of  equations  (4),  v is  conserved  in  the  mean  which  implies  that 

e = pe  + *5pu2  + k (6) 

where  k is  a constant.  The  constant  plays  no  role  and  we  are  indeed  in  the 
situation  of  system  (2).  Equations  similar  to  (4),  (5),  and  (6)  hold  in  a 
y-sweep. 

Two  questions  are  fairly  obvious  concerning  this  procedure.  Is  the  two- 
dimensional  Riemann  problem  consistently  approximated  by  the  fractional  step 
procedure?  Whether  the  answer  to  the  first  question  is  yes  or  no,  does  the 
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technique  converge  to  the  solution  of  the  general  Initial  value  problem,  assuming 
that  the  one-dimensional  Glimm  method  does  so?  Unfortunately,  the  answers  to 
these  questions  are  not  known.  The  best  that  one  can  do  is  to  obtain  numerical 
results  for  the  technique.  In  the  case  of  smooth  solutions,  a general  theory  is 
available  concerning  the  accuracy  of  operator  splitting  [Reference  12]. 

There  are  only  two  such  tests  available  at  present  in  the  literature. 

Concus  and  Proskurowski  [Reference  11]  consider  the  case  of  a single  conservation 
law;  their  initial  conditions  are  somewhat  specialized.  We  shall  not  take  up 
their  results  here  except  to  note  that  they  are  good.  The  system  (3)  has  been 
solved  by  Chorin  [Reference  3]  for  the  case  of  a channel  with  a ramp.  Figure  (7a) 
shows  the  geometry  and  the  initial  conditions  (uniform  flow).  It  is  assumed  that 
the  gas  is  polytropic.  There  are  two  adjustable  parameters,  M = Mach  number  and 
Y = ratio  of  specific  heats.  Two  cases  are  studied  for  which  an  analytical  steady 
state  solution  is  known  to  exist.  These  analytical  solutions  are  illustrated  in 
Figures  (7b)  and  (7c).  There  is  a qualitative  bifurcation  in  the  solution  from  a 
regular  to  a Mach  reflection  pattern  as  the  parameters  are  varied.  The  results  for 
this  problem  reported  by  Chorin  arise  from  using  Glimm' s method  for  the  equations  of 
unsteady  gas  dynamics  (i.e.,  system  (3))  and  waiting  for  the  approach  to  the  steady 
state.  The  results  are  very  encouraging;  the  qualitative  features  are  correct  in 
both  cases,  quantitative  results  are  very  close,  and  the  number  of  mesh  points  is 
small. 


GLIMM’ S METHOD  - EXTENSION  TO  EQUATIONS  IN  WEAK  CONSERVATION  FORM.  In  the 
case  of  one  space  dimension,  the  problem  to  be  solved  is 


u.  + f (u)  * -w(u) 

— t r 

u(r,0)  given. 


(7) 


As  before,  u represents  the  vector  of  "conserved"  quantities.  The  main  situation 
of  interest  is  the  equations  of  one-dimensional  gas  dynamics  with  spherical  or 
cylindrical  symmetry.  Explicitly,  we  have 


Pt  + (pu)r  = - (n  - 1)du/ r 

(ou)t  + (pu2  + p)r  =*  - (n  - l)pu2/r 

et  + ((e  + p)u)r  = - (n  - 1) (e  + p)u/r 
2 

e - pc  + *5pu  , c * F (p,p) 


where  r is  the  radial  coordinate,  u is  the  velocity  in  the  radial  direction,  and 
n is  a constant  which  is  equal  to  2 for  cylindrical  symmetry  and  3 for  spherical 
symmetry.  The  remainder  of  the  notation  is  as  before.  The  new  feature  here  is 
the  appearance  of  the  inhomogeneous  term  w(u);  it  is  not  possible  to  put  the 
system  (8)  in  strong  conservation  form.  Another  new  problem  is  that  w(u)  depends 
explicitly  on  r;  for  problems  involving  significant  phenomena  near  the  origin,  the 
equations  become  singular. 


12.  Gottlieb,  D. , "Strang-Type  Difference  Schemes  for  Multidimensional  Problem 
Siam  J.  Nura.  Anal.  9,  Dec.  1972,  p.  650  & 661. 
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(7a)  CONFIGURATION  FOR  REGULAR  TO  MACH  REFLECTION  TRANSITION,  a = tan*1  (1/5) 


(7b)  THE  REGULAR  REFLECTION  SOLUTION  (7c)  THE  MACH  REFLECTION  SOLUTION  FOR 


FOR  M * 2.0  and  y ■ 1.4  M - 1.6  and  y ■ 1.2 

FIGURE  7.  TRANSITION  FROM  REGULAR  TO  MACH  REFLECTION.  THE  BOTTOM  TWO  DRAWINGS  SHOW  THE 
FLOW  PATTERNS  RESULTING  FROM  ANALYTIC  SOLUTIONS  WITH  M,  y AS  INDICATED. 
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Application  of  Glimm's  method  to  the  system  (7)  involves  the  notion  of 
operator  splitting.  In  the  general  case,  given  a system  ut  = A(u)  + B (u) , one 
may  consider  solving  ut  = A(u)  and  ut  = J3(u)  in  sequence;  as  the  time  step 
approaches  zero,  one  hopes  that  the  solution  of  the  latter  problem  approaches  the 
solution  of  the  original  problem.  In  special  cases  (usually  linear  operators)  this 
can  be  proved,  for  smooth  solutions,  using  functional  analysis  (product  formulas), 
or  a simple  error  analysis  [Reference  12].  The  reason  for  using  operator  splitting 
is  that  one  can  solve  ut  = A(u)  and  ut  = B(u)  separately  but  not  the  original 
equation.  This  remark  applies  to  the  system  (7);  consider  the  splitting 


ut  + f (u)r  = 0 
u = - w(u) 


(9) 


of  (7).  The  first  half  of  (9)  is  now  in  strong  conservation  form  and  the  basic 
Glimm's  method  applies  to  its  numerical  solution.  The  second  half  of  (9)  is  a 
system  of  ordinary  differential  equations,  one  for  the  value  of  each  conserved 
quantity  at  each  mesh  point.  Observe  that  conserved  quantities  at  different  mesh 
points  are  uncoupled  in  this  system;  however,  for  the  case  of  gas  dynamics  - 
equations  (8),  the  three  conserved  quantities  at  a given  mesh  point  are  coupled. 

In  any  event,  it  is  a fairly  easy  system  to  solve  (numerically).  The  1/r  term  is 
treated  as  a constant  given  by  r = iAx  for  the  ith  mesh  point  (other  alternatives 
are  available  which  are  not  treated  here).  Note  also  that  the  origin  is  a 
boundary  condition  for  the  first  half  of  (9);  therefore,  equations  for  ii(x  = 0)  do 
not  appear  at  all  in  the  second  half  of  (9)  because  they  are  not  needed.  In  other 
words,  the  singularity  at  the  origin  has  disappeared  (at  least  numerically).  The 
ideas  involved  here  were  first  presented  by  Sod  [Reference  13]. 


The  natural  question  here  is  whether  the  splitting  (9)  converges  to  the  full 
operator  (7)  as  the  time  step  approaches  zero,  under  the  assumption  that  the 
Glimm  operator  itself  converges  (of  course,  there  is  no  question  about  the 
convergence  of  numerical  approximations  to  the  system  of  ordinary  differential 
equations).  There  are  some  positive  results  for  the  case  of  a single  conservation 
law  (private  communication  with  G.  Sod).  However,  little  is  known  for  the  case 
of  gas  dynamics,  equations  (8). 

An  alternative  to  Sod's  splitting  procedure  would  be  to  solve  the  Riemann 
problem  for  the  full  system  (7).  This  is  substantially  more  difficult  than  the 
corresponding  planar  problem,  but  it  appears  that  a mathematical  analysis  would 
yield  results — this  is  not  clear  for  the  two-dimensional  Riemann  problem  discussed 
in  the  previous  section.  Indeed,  through  the  use  of  power  series  expansions,  the 
Riemann  problem  for  the  system  (8)  of  gas  dynamics  with  spherical  symmetry 
(n  * 3)  has  been  solved  approximately  in  the  strong  shock  case  (i.e.,  the  ratio 
P^/Pj.  is  large).  Details  may  be  found  in  References  14  and  15.  Figure  8 is  a 


13.  Sod,  G.  A.,  "A  Numerical  Study  of  a Converging  Cylindrical  Shock,"  J.F.M.  83, 
1977,  p.  785-794. 

14.  Friedman,  M.  P.,  "A  Simplified  Analysis  of  Spherical  and  Cylindrical  Blast 
Waves,"  J.F.M.  11,  1962,  p.  1-15. 

15.  Holt,  M. , "The  Initial  Behavior  of  a Spherical  Explosion.  I.  Theoretical 
Analysis,"  Proc.  of  the  Royal  Society,  A,  234,  1956,  p.  89-109. 
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schematic  illustration  of  the  solution  in  the  strong  shock  (moving  away  from  the 
origin)  case.  One  sees  immediately  that  the  solution  is  no  longer  simply  a 
function  of  r/t,  but  depends  on  r,t  separately.  This  makes  the  application  of 
this  solution  to  Glimm's  method  awkward  and  inefficient.  Also,  for  such  an 
application  to  be  considered,  it  will  be  necessary  to  extend  the  analysis  to  the 
entire  spectrum  of  pressure  ratios;  however,  pjt/pr>>l  is  probably  the  most 
difficult  case.  Hence,  for  the  time  being,  the  splitting  algorithm  is  the  only 
way  to  apply  Glimm's  method  to  systems  in  weak  conservation  form.  Even  if  tech- 
niques involving  exact  or  approximate  solutions  of  the  full  Riemann  problem  (e.g.. 
Figure  8)  prove  infeasible,  it  would  be  useful  to  have  these  solutions  available 
as  a test  problem  for  the  splitting  method. 

We  now  present  the  available  numerical  evidence  for  the  splitting  algorithms 
applied  to  the  system  (8),  one-dimensional  gas  dynamics.  Sod  studies  [Reference 
13]  the  case  of  a converging  cylindrical  shock  induced  by  a cylindrical  shock 
tube  problem;  this  problem  has  an  extensive  literature  from  both  the  computational 
and  experimental  sides.  The  singularity  at  the  origin  (pressure  approaches 
infinity  as  the  shock  approaches  the  origin)  makes  this  a difficult  problem  for 
standard  finite  difference  schemes.  The  results  presented  by  Sod  are  in  good 
agreement  with  experimental  data  and  are  an  improvement  on  prior  computations. 

The  method  handles  the  singularity  accurately  and  without  difficulty,  the  correct 
sequence  of  waves  is  formed;  wave  interactions  appear  correctly,  and  the  reflection 
of  waves  from  the  origin  is  as  expected.  In  short,  the  results  for  this  problem 
have  all  of  the  desirable  properties  that  we  have  seen  for  Glimm's  method  in  the 
plane  case,  despite  the  added  complication  of  the  nontrivial  symmetry  terms. 

In  the  case  of  spherical  symmetry,  the  Primakoff  problem  provides  an  exact 
solution.  The  derivation  of  this  similarity  solution  may  be  found  in  Reference  7. 
The  solution  consists  of  a curved,  nonuniform  waveform,  which  decays  to  zero  at 
the  origin,  behind  an  expanding  shock  wave  front  which  is  sharply  peaked. 
Analytically,  the  solution  is  given  by 

u(r,t)  = O.lrt  p(r,t)  = pQr\  12^/25k;  p(r,t)  = 4pQrt  2^5/3k  (10) 

where  pQ  is  the  constant  ambient  density  and  k is  a constant.  The  location  of  the 
shock  is  given  by 

R(t ) = kt2/5.  (11) 

We  chose  (arbitrarily)  t = 1 for  the  initial  time  in  our  computations  and  used  the 
discretized  versions  of  equations  (10)  and  (11)  for  initial  conditions.  At 
t * 5.6,  R(t)  is  roughly  double  R(l)  and  so  this  was  chosen  as  the  final  integra- 
tion time.  Figures  9 and  10  reproduce  density  and  pressure  profiles,  respectively, 
at  t = 5.6,  for  both  the  exact  and  numerical  solution.  Figure  11  is  the  peak 
pressure  vs.  distance  profile  for  the  wave  from  t = 1 to  t = 5.6;  this  is  an 
important  function  in  evaluating  blast  wave  effects.  It  can  be  seen  that  the 
numerical  results  for  Glimm's  method  are  in  agreement  with  the  exact  solution;  in 
the  run  ^presented  here  the  shock  location  is  exact,  but  in  other  runs  it  is  off  a 
few  mesh  points.  It  is  not  clear  at  this  point  whether  this  error  is  caused 
by  the  random  choice  element  or  by  the  operator  splitting;  in  any  event, 
we  can  report  that  reducing  the  mesh  will  reduce  the  error  in  shock 
location  indicating  that  convergence  is  obtained  for  this  problem.  Using  a new 
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idea  reported  in  Glimm  & Marchesin^  (1978),  it  may  be  possible  to  track  the  shock 
front  (almost)  exactly,  yet  within  the  context  of  the  random  choice  method;  we 
have  not  tested  these  techniques  as  yet. 


The  case  of  two  dimensional  axisymmetric  compressible  flows  can  be  handled  by 
combining  the  fractional  step  approach  of  the  last  section  with  the  operator 
splitting  of  this  section.  The  equations  are  of  the  form 

+ f(u)r  + &(u)2  = ~w(u).  (12) 

The  idea  is  to  use  operator  splitting  and  solve  in  sequence 


uc  + f(u)r  + £(u)z  = 0 
u = - w(u) . 


(13) 


The  first  equation  is  solved  using  the  fractional  step  approach  of  Chorin; 
the  second  equation  is  simply  a system  of  ordinary  differential  equations.  At 
the  present  time,  there  is  no  published  work  available  on  Glimm's  method  applied 
to  this  situation. 


Operator  splitting  can  also  be  applied  to  the  equations  of  gas  dynamics  with 
chemical  reaction.  In  this  case,  the  vector  w(u)  in  equation  (7)  or  (12)  will  be 
zero  except  in  the  last  component  where  a term  for  the  net  energy  gain  from 
chemical  reaction  will  appear.  Of  course,  it  would  be  possible  to  combine 
chemical  reaction  with  nontrivial  spatial  symmetries  by  including  a more  complicated 
vector  w(u) . The  problem  with  the  operator  splitting  approach  is  that  the  time 
scales  for  hydrodynamics  and  chemistry  are  (typically)  off  bv  several  orders  of 
magnitude  leaving  the  investigator  with  the  choice  of  either  using  an  extremely 
small  time  step,  or  using  a very  inaccurate  integration  of  the  chemical  phenomena. 
(This  dilemma  is  inherent  in  the  problem,  not  just  the  Glinm  approach  to  the 
problem.)  An  ingenious  technique  avoiding  this  difficulty  is  presented  by  Chorin 
[Reference  6];  the  idea  is  to  explicitly  solve  Riemann  problems  with  chemistry 
using  the  assumption  that  any  reaction  that  takes  place  during  a time  step  takes 
place  instantaneously  at  the  beginning  of  the  time  step.  This  preserves  the  self- 
similar nature  of  the  solution  to  the  Riemann  problem,  which  now  may  consist  of 
Chapman-Jouguet detonat ions  and/or  strong  detonations  as  well  as  hydrodynamic  shocks, 
rarefaction  waves,  and  contact  discontinuities.  (The  deflagration  case  can  be 
handled  as  well,  but  a heat  conduction  term  must  be  added  to  the  system  of  equa- 
tions on  physical  grounds.)  In  the  paper,  numerical  evidence  is  presented  for 
both  strong  and  Chapman-Jouguet  detonations  with  excellent  results.  However,  the 
reaction  rate  used  is  in  close  agreement  with  the  unphysical  hypothesis  necessary 
in  Chorin's  construction  of  the  solution.  Thus,  further  work  is  necessary  in  this 
area,  but  it  does  appear  that  Glimm's  method  can  be  used  rather  efficiently  in 
problems  involving  chemical  reaction. 

APPLICATIONS  TO  BLAST  WAVE  PHENOMENA.  In  the  field  of  numerical  modeling  of 
blast  waves,  one  may  perceive  two  long-range  goals.  The  first  is  to  develop 
methods,  based  on  the  governing  partial  differential  equations,  which  are 
sufficiently  accurate  and  reliable  that  experimental  confirmation  of  the  results 


16.  Glimm,  J.  and  Marchesin,  D.,  "A  Random  Numerical  Scheme  for  One  Dimensional 
Fluid  Flow  with  High  Order  of  Accuracy,"  preprint.  May,  1978 
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Is  not  required;  the  achievement  of  this  goal  has  obvious  practical  importance. 

The  second  goal  is  to  develop  the  methods  to  such  an  extent  that  they  are  useful 
in  a theoretical  analysis  of  blast  wave  phenomena — a field  of  study  which  is  by 
no  means  completely  understood.  These  two  goals  are  clearly  interrelated; 
however,  the  exact  relationship  is  not  obvious  at  this  time. 

The  difficulties  in  solving  blast  wave  problems  can  usually  be  traced  to  the 
existence  of  discontinuities  in  the  flow  field — this  includes  shock  waves,  slip 
surfaces,  and  material  interfaces.  The  calculation  in  a neighborhood  of  the 
discontinuities  is  both  the  most  difficult  part  to  solve  numerically  and  the  most 
important  part  from  a practical  point  of  view.  Indeed,  the  solution  is  almost 
useless  unless  it  contains  peak  pressure  vs.  distance  plots,  shock  trajectories, 
interface  trajectories,  etc.  Furthermore,  a typical  phenomenon  in  blast  wave 
computations  is  the  formation,  interaction,  and  disappearance  of  discontinuities  in 
the  midst  of  the  computation;  this  includes  bifurcations  of  discontinuity  patterns 
such  as  the  transition  from  regular  to  Mach  reflection.  The  physical  effects  of 
this  behavior  are  critical  in  the  applications. 

The  development  of  numerical  techniques  to  handle  these  problems  began  in 
World  War  II  and  continues  today.  Reference  17  contains  a recent  review  of 
several  of  these  methods.  Most  rely  on  artificial  viscosity  and  "shock  capturing" 
while  some  jse  shock  fitting  either  with  or  without  the  usual  techniques  in  other 
parts  of  the  flow  field.  The  basic  problem  with  the  former  method  is  that  it 
smears  out  shocks  over  several  mesh  points;  the  resulting  distortion  near  the 
front  is  very  undesirable  in  itself  and  can  easily  lead  to  further  errors  in 
interaction  problems.  On  the  other  hand,  shock  fitting  would  be  ideal  for  the 
problems  at  hand  if  it  could  be  implemented  in  a straightforward  fashion.  However, 
the  necessary  theory  is  not  available  to  track  multidimensional  discontinuities 
(in  general,  curved)  and  their  interactions  and  transitions.  Even  if  this  were 
feasible,  it  would  be  still  necessary  to  anticipate  discontinuity  formation 
(either  automatically  or  otherwise)  and  prepare  special  algorithms  for  each  along 
with  the  interactions  which  occur;  this  can  be  a very  ad  hoc  procedure. 

The  advantages  of  Glimm's  method  are  clear;  discontinuities  are  tracked 
without  smearing,  and  the  formation,  interaction,  and  transition  of  discontinuities 
is  computed  automatically.  This,  at  least,  is  the  situation  in  one  space 
dimension.  The  most  important  open  question  for  the  method  is  its  utility  in 
higher  dimensional  problems  (this,  of  course,  is  also  the  case  for  most  other 
numerical  techniques  in  compressible  fluid  flow).  The  results  of  Chorin  presented 
in  the  section  on  the  fractional  step  approach  are  extremely  promising  in  this 
regard;  they  indicate  that  the  sharp  resolution  of  discontinuities  of  the  one- 
dimensional scheme  is  also  a property  of  the  two-dimensional  fractional  step 
approach. 

We  intend  to  apply  the  method  to  a variety  of  multidimensional  blast  wave 
situations  beginning  with  the  problem  of  a spherical  or  cylindrical  charge  at 
height  of  burst.  This  becomes  axisymmetric  (two-dimensional)  at  the  moment  of 
impact  of  the  shock  on  the  ground.  This  problem  contains  transitions  from  regular 
to  Mach  reflection  in  addition  to  the  initial  strong  shock  wave  followed  by  a 
contact  discontinuity.  If  the  initial  energy  is  high  enough,  there  is  a further 
transition  to  multiple  Mach  reflection — this  is  amply  documented  by  experiment 
although  not  understood  theoretically.  In  a closely  related  problem,  we  will 
study  the  two-dimensional  plane  ramp  problem  with  an  incident  plane  shock  as 

17.  Sod,  0.  A.,  "A  Survey  of  Several  Finite  Difference  Methods  for  Systems  of 

Nonlinear  Hyperbolic  Conservation  Laws,"  J.C.P.  27,  April  1978,  p.  1-31. 
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initial  data  (this  is  the  "shock  tube"  problem  whereas  Chorin  studied  the  "wind 
tunnel"  problem),  see  Figure  7.  There  is  a large  body  of  experimental  data 
available  for  this  situation  and  the  resulting  unsteady  flow  fields  are  similar 
to  the  axisymmetric  charge  at  height  of  burst  problem.  For  the  final  applications, 
it  will  be  necessary  to  use  realistic  equations  of  state;  the  appendix  indicates 
that  this  will  be  feasible.  This  will  also  be  very  important  in  the  more  difficult 
charge  at  height  of  burst  above  water  problem.  In  addition  to  the  physical 
phenomena  listed  above,  this  problem  also  contains  refraction  of  discontinuities 
as  well.  An  advantage  of  Glimm's  method  for  this  problem  is  that  as  the  leading 
shock  approaches  the  water,  there  is  no  reason  for  the  time  step  to  approach  zero 
as  is  the  case  for  any  Lagrangian  scheme  used  on  the  problem;  as  the  ambient  atmos- 
phere gets  squeezed  between  the  shocked  air  and  the  water,  the  probability  of  a ran- 
dom number  lying  in  that  portion  of  the  Riemann  problem  solution  corresponding  to 
ambient  atmosphere  will  approach  zero  and  eventually  the  ambient  atmosphere  must 
disappear.  Since  water  is  only  slightly  compressible,  its  equation  of  state  will 
not  even  resemble  that  for  a y-law  gas  and  the  development  of  a Riemann  problem 
solver  for  water  will  be  important. 

Other  problems  involving  material  interfaces  include  a charge  impacting  on  a 
solid  surface  and  the  design  of  shaped  charges  (which  include  a metal  liner  inside 
the  charge).  On  the  one  hand,  these  problems  should  be  amenable  to  the  approach 
embodied  in  Glimm's  method  and  the  advantages  of  the  latter  indicated  above  will 
be  important.  On  the  other  hand,  the  method  has  not  yet  been  tested  in  problems 
of  this  type  involving  nonreflecting  material  interfaces.  In  one  space  dimension, 
this  type  of  complication  is  not  a problem;  indeed,  Chorin  [Reference  3]  points 
out  that  interfaces  can  be  tracked  with  the  same  accuracy  (in  the  mean)  as  can  be 
attained  by  the  method  for  shocks.  In  higher  dimensions,  this  property  remains  to 
be  tested. 

We  note  also  that  research  is  ongoing  to  apply  Glimm's  method  to  problems 
involving  combustion  and  detonation.  We  have  outlined  this  development  in  the 
preceding  section.  Thus,  it  may  be  possible  in  the  near  future  to  use  the  method 
to  model  the  detonation  wave  from  an  explosive  instead  of  simply  using  experimental 
data  to  set  up  initial  conditions  for  a hydrodynamic  shock  wave. 

Finally,  it  is  possible  to  couple  Glimm's  method  to  other  techniques  in 
separate  flow  regions.  Chorin  studies  the  interaction  of  a laminar  incompressible 
boundary  layer  with  an  inviscid  compressible  interior  flow  [Reference  18].  Of 
course,  a completely  different  method  is  necessary  for  the  boundary  layer.  Glimm 
and  Marchesin  [Reference  16]  consider  "random  shock  tracking."  This  adaptability 
of  the  method  may  prove  useful  in  future  applications. 


18.  Chorin,  A.  J. , "Vortex  Sheet  Approximation  of  Boundary  Layers,"  J.C.P.  27, 
June  1978,  p.  428-442 
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APPENDIX 

THE  GODUNOV  ITERATION  FOR  SOLVING  RIEMANN  PROBLEMS 


Figure  1 contains  most  of  the  notation  which  we  shall  use.  In  addition,  set 
t - 1/p,  e = internal  energy,  and  U = shock  speed.  Solving  a Riemann  problem 
means:  given  u^  = (pj.u^.p^)  and  ur  = (pr,ur,pr),  calculate  P*.u*,p*£,P*r,  the 
slopes  of  the  shock  (or  shocks),  contact  discontinuity,  and  the  leading  and 
trailing  edges  of  the  rarefaction  fan  (or  fans),  and  the  solution  inside  the  fan. 
In  the  application  to  Glimm's  method,  it  is  never  necessary  to  calculate  the 
entire  solution;  one  always  finds  p*  and  u*  and  certain  other  parts  of  the 
solution  depending  on  the  random  number  selected.  In  particular,  it  is  never 
necessary  to  explicitly  find  the  entire  rarefaction  fan. 

The  first  step  in  the  solution  is  the  Godunov  iteration  for  the  determination 
of  p*.  This  will  be  presented  in  some  detail  for  a general  equation  of  state. 
Next,  the  explicit  formulas  in  the  iteration  will  be  derived  for  a y-law  gas  and 
the  complete  solution  of  the  Riemann  problem  will  be  set  down  in  this  case.  This 
material  is  taken  from  References  3,  6,  7,  18,  and  19.  Finally,  the  Godunov 
iteration  will  be  repeated  for  more  complicated  equations  of  state  which  arise  in 
the  study  of  chemical  explosives.  It  will  follow  from  this  analysis  that  Glimm's 
method  can  be  extended  in  a straightforward  fashion  to  problems  involving  general 
equations  of  state.  The  remaining  question  is  the  (numerical)  efficiency  of  the 
method  for  complicated  equations  of  state.  This  is  not  resolved  here  except  for 
a few  general  comments.  In  particular,  we  do  not  address  the  possibility  of 
replacing  Godunov's  iteration  scheme  (which  is  very  efficient  in  the  y-lav  gas 
case)  by  some  other  scheme.  It  will  be  clear  that  the  iteration  for  p*  is  the 
bulk  of  the  computing  time  in  the  application  of  Riemann  problem  solutions  to 
Glimm's  method.  Thus,  we  will  not  repeat  the  derivation  of  the  complete  solution 
to  the  Riemann  problem  for  these  equations  of  state. 

Given  p*  and  pr,  it  follows  that  the  *-state  and  the  right  state  are 
separated  by  either  a shock  wave  or  a centered  rarefaction  fan  according  as 
p*  * pr  or  p*  * pr,  respectively  (entropy  condition).  The  same  result  holds  for 
pA  and  p^.  So,  in  order  to  proceed  in  the  iteration,  we  set  down  the  governing 
equations  for  a shock  wave  and  for  an  isentropic  rarefaction  wave.  To  begin, 
suppose  that  the  equation  of  state  is  of  the  form 


19.  Sod,  G.  A.,  "The  Computer  Implementation  of  Glimm's  Method,"  Lawrence 
Livermore  Laboratory  Report  UCID-17252,  1976 
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P " p(p,e),  (Al) 

and  that  equation  (Al)  may  be  solved  explicitly  for  the  internal  energy  e, 

e = e(p,p).  (A2) 

By  using  the  thermodynamic  identity 

(|f)s  - CS9,  - p(|f)T  <«> 

and  the  equation  of  state  (Al)  and  (A2),  the  equation  for  an  isentrope  can  always 
be  put  in  the  form 

^ = h (t , p)  (A4) 

or 

p = p(t)  - /T  h (t ' , p ' )dr ' . (A5) 


The  solution  of  (A5)  may  or  may  not  be  obtainable  in  closed  form;  in  any  event, 
this  solution  will  contain  a constant  of  integration  which  depends  on  the  isentrope. 
Additionally,  there  are  also  Riemann  invariants  which  are  of  the  form 


where 


Up)  ± u = const. 


(A6) 


Up) 


c 'd 


= / 


P’  d^_ 


P c' 


(A7) 


The  constants  of  integration  in  (A7)  may  always  be  taken  to  be  zero  for  gases; 
see  Reference  7.  The  constant  in  (A6)  depends  on  the  value  of  the  entropy.  The 
sound  speed  c may  be  determined  from 


2 

c 


(A8) 


the  right-hand  side  of  equation  (A8)  is  obtainable  from  equation  (A3).  Therefore, 
for  the  case  of  a rarefaction  wave  between  the  *-state  and  the  right  state, 

*(Pr)  - ur  = UP*)  - u*  (A9) 

H(Pr,Pr)  = H(p*,p*)  (A10) 

where  H is  obtained  from  integration  of  equation  (A5).  Equations  (A9)  and  (A10) 
follow  because  the  rarefaction  wave  is  isentropic.  Now,  consider  the  case  of  a 
shock  wave  between  the  *-state  and  the  right.  We  have  the  Rankine-Hugoniot  condi- 
tions which  are  completely  independent  of  the  equation  of  state: 

Pr(ur  - U)  = P*(u*  - U)  « - M (All) 

0r(ur  " U>2  + Pr  - P*(u*  “ u)2  + P*  (A12) 
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er  - e*  - ((r*  - Tr)/2)(p*  + pr)  (A13) 

From  equations  (All)  and  (A12),  It  follows  that 

M2  = - (p*  - Pr)/(T*  - Tr).  (A14) 

The  exact  analogues  of  equations  (A9)  - (A14)  hold  for  the  waves  between  the  left 
state  and  the  (*)-state. 

After  an  initial  guess  for  p*  has  been  made,  the  Godunov  iteration  proceeds 
by  defining 


Mr  = (pr  - p*)/(ur  - u*)  (A15) 

M,,  “ - (Pe  - P*)/(ua  - u*).  (A16) 

Observe  that  equations  (A15),  (Alb)  imply 

P*  - K ~ ur  + Pr/Mr  + PJl/Me)/((l/Mr)  + (1/M^)).  (A17) 

Given  an  iterate  p*,  the  Godunov  iteration  updates  Mr  and  M to  obtain  the 
next  iterate  p*"*^.  This  is  repeated  until  convergence  is  reached.  Convergence 
criteria,  accelerated  convergence  techniques,  and  the  optimal  choice  of  the  initial 
guess  pO  will  not  be  discussed  here;  see  Reference  6.  The  quantities  Mr  and  Mj^ 
cannot  be  used  in  the  form  given  in  (A15),  (Alb)  for  this  process  because  they 
involve  the  unknown  u*.  The  main  idea  in  the  iteration  is  to  rewrite  Mr  and  as 
functions  of  (pr,p  ,p*),  (P£*P{.,P*)  respectively,  using  equations  (Al)  - (A14) . 

This  is  now  carried  out  for  the  wave  connecting  the  *-state  with  the  right  state. 

The  wave  moving  to  the  left  may  be  handled  similarly. 

Suppose  that  this  wave  is  a shock.  From  equation  (All),  ur  - u*  = - M(xr  - t*) 
which  implies  that  Mr  = M by  equation  (A14) . This  means  that 

M2  = - (pr  - P*)/(Tr  - t*).  (A18) 

This  eliminates  u*  but  introduces  t*  into  the  equation  for  Mr.  However,  we  can 
obtain  t*  ■ T*(pr»Pr»P*)  by  using  the  Rankine-Hugoniot  condition  (A13)  after 
expressing  er  - er(pr,pr),  e * = E*(P*.P*)  through  the  equation  of  state  (A2).  This 
involves  the  solution  of  a nonlinear  algebraic  equation.  We  are  assuming  that  the 
equation  of  state  is  such  that  this  nonlinear  equation  has  exactly  one  root  which 
satisfies  the  governing  physical  laws.  Inserting  this  expression  into  equation 
(A18) , we  have  Mr  * Mr(pr,p*,pr)  as  desired. 

Now  suppose  that  the  right  wave  is  an  isentropic  rarefaction.  Then,  equations 
(A15)  and  (A9)  lead  directly  to 

Mr  - (Pr  - P*)/a(Pr)  - Up*))  (A19) 

It  follows  from  equations  (A2),  (A7),  and  (A8)  that  fc(pr)  is  a function  of  pr  and 
pr,  and  that  Up*)  is  a function  of  p*  and  p*.  The  remaining  task  is  to  eliminate 
P*.  This  may  be  accomplished  by  solving  the  isentropic  law  (A10)  for  p*  in  terms 
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°f  pr,pr,  and  p*.  Once  again,  this  is  a nonlinear  algebraic  equation  and  we  make 
the  same  assumption  about  its  solution  as  was  made  above.  Having  solved  this 
equation,  one  sees  that  equation  (A19)  exhibits  the  quantity  Mj.  in  the  appropriate 
form. 


Specializing  to  the  case  of  a y-law  gas,  the  equation  of  state  is 


1 £ 

c = — - . 

Y-l  P 

Application  of  the  thermodynamic  identity  (A3)  leads  to 

dp/dr  = -yp/t 

and  the  solution  of  this  ordinary  differential  equation  is 

PP  y = A(S) 


(A20) 


(A21) 


(A22) 


where  A(S)  is  a constant  of  integration  (which  depends  on  the  entropy,  S).  This  is 
the  isentropic  law,  and  equation  (A10)  takes  the  form 


P P 
*r  r 


P*P 


-Y 


(A23) 


Similarly,  one  uses  equations  (A3),  (A8),  and  (A20)  to  show  that  the  sound  speed 
satisfies 


YP/P- 


(A24) 


By  inserting  equation  (A24)  into  equation  (A7)  for  £(p)  and  using  equation  (A22), 
we  see  that  equation  (A9) , which  expresses  the  constancy  of  the  right  Riemann 
invariant  across  the  right  wave,  may  be  put  in  the  form 


2c* 


2c 


u = ; U . . 

Y-l  r Y-l  * 


(A25) 


Therefore,  if  the  right  wave  is  a rarefaction,  equations  (A15) , (A24) , and  (A25) 
lead  to 


",  ■ ^ [£)”  - £>”] 


(A26) 


Using  the  isentropic  law  (A23),  it  follows  (after  some  algebra)  that 

y-l/2y\ 


(A27) 


Thus,  the  general  case  (A19)  for  Mj.  has  reduced  to  a simple  algebraic  formula. 
Suppose  now  that  the  right  wave  is  a shock.  Using  the  equation  of  state  (A20)  and 
the  Ranklne-Hugoniot  condition  (A13),  one  obtains 
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where  u 


2 


Y-l 

Y+1‘ 


T, 


* 


Wp,\ 

p2pr  * p*/ 


When  substituted  Into  equation  (A18), 


M 

r 


(Vr) 


+ 


I±1 

2 


(A29) 


Equation  (A29)  is  of  the  required  form.  Once  again,  it  is  merely  an  algebraic 
expression  for  Mr.  Collecting  our  results,  we  see  that 


where 


Mr  “ (PrPrK*(P*/Pr) 


(A30) 


<K  a) 


Izl  ..  _ < , 

2,’i  - 


(A31) 


Observe  that  is  continuous  at  a = 1;  indeed  C 1 ) = Y‘.  Formula  (A30)  also  holds 
for  Mjj  except  that  pr  is  replaced  everywhere  by  p?.  Thus,  formula  (A30) , (A31), 
and  (A17)  define  the  iteration  for  p*  in  the  case  of  a y-law  gas.  The  question  of 
convergence  will  not  be  taken  up  here  except  to  note  that  the  iteration  converges 
very  quickly  in  practice.  Reference  6 may  be  consulted  for  further  details. 


The  remainder  of  the  solution  to  the  Riemann  problem  is  now  derived.  The 
equations  will  be  for  a y-law  gas,  but  the  procedure  easily  generalizes  to  an 
arbitrary  equation  of  state.  One  notes  immediately  that 


u*  " (P*  " Pr  + Mrur  + V*)/(Mr  + V (A32) 

upon  elimination  of  p^  from  (A15) , (A16) . Observe  that  the  equation  for  the  slip 
line  is  dx/dt  - u*.  Suppose  there  is  a shock  wave  on  the  right  (i.e.,  p*  > pr). 
Equations  (All)  gives  us  the  shock  speed  U and  the  density  to  the  riRht  of  the 
slip  line  p*r  since  the  quantity  M is  known  from  the  iteration  for  p*.  A left- 
facing shock  is  handled  similarly.  The  remaining  case  is  a rarefaction.  Suppose 
it  is  the  right  wave.  The  rarefaction  is  then  bounded  on  the  right  by  the  line 
dx/dt  ■ u_  + cr  - u_  + (yppTr)  and  so  this  line  is  known  a priori.  It  is  bounded 
on  the  left  by  the  line  dx/dt  - u*  + c*  - u*  + (yp*t*)H.  The  quantity  x*  (i.e., 
l/p*r)  is  not  known  but  can  be  found  immediately  from  equation  (A25).  Thus,  the 
left  slope  and  the  density  P*r  are  solved  for  simultaneously.  If  the  solution 
is  desired  at  some  point  p inside  the  rarefaction,  one  equates  the  slope  of  the 
characteristic  dx/dt  * u+c  to  the  slope  of  the  line  connecting  p and  the  origin. 
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This  is  correct  because  the  rarefaction  fan  is  centered  at  the  origin.  Hence, 
u+c  = 6 where  8 is  the  slope.  Also,  the  constancy  of  the  right  Riemann 
invariant  and  the  isentropic  law  are  available  to  give  two  more  equations  for 
the  unknowns  p,p,u.  Hence,  there  is  a total  of  three  equations  in  three  unknowns 
and  the  solution  can  be  found.  Observe  that  knowledge  of  the  *-state  is  not 
necessary  to  find  the  solution  inside  the  fan  after  one  knows  what  the  fan 
boundaries  are  (which  does  require  information  from  the  *-state). 


In  summary,  it  has  been  shown  that  the  general  solution  of  the  Riemann 
problem  can  be  reduced  to  simple  algebraic  formulas  in  the  case  of  a y-law  gas. 

For  applications  to  blast  wave  problems,  there  will  be  many  situations  when  a 
y-law  approximation  to  the  equation  of  state  will  be  inappropriate  either  for  the 
detonation  products,  the  ambient  atmosphere  (e.g.,  after  radiant  heating  has 
raised  the  temperature  to  a very  high  level),  or  some  other  media  in  the  problem 
(e.g.,  in  a blast  above  water,  the  water  should  be  treated  as  a compressible 
medium  if  the  blast  strength  is  great  enough).  Hence,  it  is  important  to  be  able 
to  efficiently  solve  Riemann  problems  for  more  complex  equations  of  state  in  order 
to  be  able  to  apply  Glimm's  method  to  such  problems.  In  addition,  it  is  necessary 
to  be  able  to  solve  Riemann  problems  across  a material  interface  (e.g.,  the  detona- 
tion product-ambient  atmosphere  interface).  In  this  situation,  the  left  and 
right  states  will  have  different  equations  of  state. 

Ritter  [Reference  20]  considers  an  explosive-metal  system.  For  the 
compressible  metal  plate,  Murnaghan's  equation  of  state  is  used — 


e = (v-1)  1[ (p+vb)T-vbTQ] 


(A33) 


where  v > 1,  b,  tq  are  constants.  The  isentropic  law  is  given  by  p = b[(— ) -1] 

2 P° 
and  c - v(p+b)t  gives  the  sound  speed.  An  analysis  similar  to  that  leading  up 

to  equation  (A36)  shows  that  for  the  equation  of  state  (A33),  the  quantity  Mr  is 

given  by 


Mr  - [(Pr  + b)pr]  ><K(p*+b)/(pr  + b)) 


where 


(A34) 


$(a)  = 


/ V+~ 1 V—  A.  ' 

(~T“  a + 


v-1 

2v 


1-a 

1_a(v-!)/2v 


a > 1 

a < 1 


(A35) 


Hence,  this  equation  of  state  is  no  more  difficult  to  work  with  than  that  of  a 
y-law  gas.  Glimm's  method  for  this  problem  is  compared  with  previous  experimental 
results  in  Reference  20. 


20.  Ritter,  Z.  W. , "A  New  Method  for  Calculating  Hydrodynamic  Behavior  of  Plane 
One-Dimensional  Explosive-Metal  Systems,"  1977,  preprint 
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Kury,  et  al  [Reference  21]  use  the  equation  of  state 


-R-i  t -R0x 

p A<1  - iy)e  + 8(1  - R^,e  + 7" 


(A36) 


for  a chemical  explosive  known  as  Comp  B,  Grade  A.  In  equation  (A36) , the 
quantities  A,  B,  R^,  R2*  w are  constants.  Applying  equation  (A3),  the  following 
ordinary  differential  equation  is  obtained  for  the  isentrope  (after  some  algebraic 
manipulations  using  the  equation  of  state  (A36)): 


dp  l+u  t ,arH  „ . 1T  , n/W+l  „ 2T 

df  + ~ p = A(_r " Rl)e  + B(~  " Ve 

The  integrating  factor  for  equation  (A37)  is  x^1^^  and  one  obtains  after  a 
straightforward  calculation  that  the  solution  of  (A37)  is 

p = Ae_RlT  + Be""2'  + C r<“+1). 


(A37) 


(A38) 


Here,  C is  a constant  of  integration  which  plays  the  same  role  as  the  constant 
A(S)  in  equation  (A22)  for  the  y-law  gas  isentropes;  that  is,  C is  different  for 
different  isentropes.  Combining  equations  (A8)  and  (A37),  it  is  easy  to  see  that 
the  sound  speed  is  given  by 


- -R.x  -R2t 

c = (l+w)px  - Ar(l+u  - R^x)e  - Bx(l+w  - R2x)e  ; (A39) 

substituting  the  expression  (A38)  for  the  pressure  p in  equation  (A39) , one 
obtains 

c2  = Dx’“  + AR1T2e  1 + BR2t2e  2 (A40) 


where  D is  another  constant  depending  on  the  isentrope.  Equation  (A7)  for  the 
quantity  R(o)  appearing  in  the  definition  of  Riemann  invariant  becomes 

-I  r\  _!  R,/P  R«/Pl 

i(p’)  = f°  (Dp  Z + ARjP  e"  + BR2p'^e  )2dp.  < 


(A41) 


(Note:  the  constant  w is  always  greater  than  zero  so  that  J7  (p ) does  not  blow  up 

near  p » 0.  Indeed,  according  to  Courant  and  Friedrichs  [Reference  7],  one  may 
take  8. (p ) = 0 for  p = 0,  at  least  for  gases.)  The  analogue  of  equation  (A10)  for 
the  isentropic  law  (A38)  is 

-R,  r -R0t  ,,  -R.t*  -R.t* 

w+1,  . 1 r _ 2 r.  w+1,  1 * „ 2 * 

t (p  - Ae  - Be  ) = x (pA  - Ae  - Be  ).  (A42) 

r r * x 


21.  Kury,  J.  W. , et.al.,  "Metal  Acceleration  by  Chemical  Explosions,"  Fourth 
Symposium  (Int.)  on  Detonation,  Oct.  1965,  p.  3-13. 
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Hence,  If  the  right  wave  Is  a centered  rarefaction,  Mr  Is  found  from  equation  (A19) 
which,  in  view  of  equation  (A41) , takes  the  form 


= (Pr  " P*)/ 


/ 


P*(Pr.Pr,P*) 


(Dp 


co— 2 


AR^p 


-4  “Rl/p 
e 


br2p 


•4  "R2/p  V, 

e ) dp 


(A43) 


where  p * = P*(pr.Pr»P*)  is  the  solution  of  equation  (M2).  On  the  other  hand,  if 
the  right  wave  is  a shock,  Mr  is  found  from  equation  (A18) , 


= (Pr  " P*)/(Tr  ~ T*(Pr>pr>P*)) » (A44) 

where  x*  is  the  solution  of  the  Rankine-Hugoniot  condition  (A13);  the  quantities 
Er,£*  are  easily  eliminated  in  favor  of  pr,Pr  and  p*,p*  by  direct  use  of  the 
equation  of  state  (A36)  and  the  result  is  a nonlinear  algebraic  equation  for  x* 
in  terms  of  pr,Pr,  and  p*.  Similar  results  hold  for  the  left  wave.  In  summary, 
solving  the  Riemann  problem  for  Comp  B,  Grade  A will  involve  the  solution  of  a 
nonlinear  algebraic  equation,  either  (A42)  or  (A13),  and  possibly  a quadrature 
(equation  A(43)),  at  each  step  in  the  Godunov  iteration. 


Thus,  Comp  B,  Grade  A is  an  example  of  an  equation  of  state  for  which 
nontrivial  numerical  procedures  must  be  introduced  at  each  step  in  the  iteration. 
However,  it  does  not  represent  the  worst  possible  case — in  general,  the 
differential  equation  for  the  isentrope  (A4)  will  not  be  reducible  to  quadratures 
and  it  will  be  necessary  to  use  a numerical  ordinary  differential  equation  solver 
in  order  to  obtain  the  isentropic  law.  The  situation  for  the  equation  of  state 
of  water  illustrates  this  eventuality.  Among  other  possibilities,  we  will 
consider  the  y-law  type  and  the  Sternberg-Walker  equations  of  state  which  are 
both  presented  in  Reference  22.  The  former  is  given  by 


[Pqt0  + yA(xp-T ) + (y-1)e]  (y-1)e/t  (A45) 

where  Y»A,po»to  are  constants.  A computation  shows  that  the  isentropic  law  (A4) 
takes  the  form 

^ + [(y-IHPqTq  + yA(Tq-t))  + l]p/x 

= %0/t2  ~ Ya/ t ) [~ (Pqt0  + yA(Tq-t))  ± V(p oTo+YA(T0_T>)2+4pT  ^ ’ (A46) 

The  square  root  arises  because  the  internal  energy,  e,  must  be  eliminated  using 
the  equation  of  state  (A45)  which  exhibits  an  e2  term  in  the  equation  for  p. 

Thus,  in  an  iteration  of  the  Godunov  method  which  contains  a rarefaction  as  the 
right  wave,  the  only  way  to  find  the  value  of  x*  corresponding  to  Pr,rr  and  the 
updated  p^  is  to  numerically  solve  the  ordinary  differential  equation  (A46)  using 


22.  Enig,  J.  W. , "The  Unsteady  Regular  and  Mach  Reflection  Resulting  from  the 
Interaction  of  Spherical  Explosion  Shock  Waves  in  Water,"  Sixth  Symposium 
(Int.)  on  Detonation  1976,  p.  570-590. 
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P = pr  at  t = Tr  as  the  Initial  data  until  p = p*  is  reached  at  which  point  the 
integration  is  stopped  and  x*  is  set  equal  to  the  current  value  of  x.  Also,  one 
can  proceed  as  in  the  previous  example  and  obtain  the  sound  speed  as  a function  of 
p and  x by  combining  equations  (A8)  and  (A46) ; unfortunately,  equation  (A46)  must 
be  used  to  eliminate  p in  terms  of  r which  means  that  the  equation  for  the  Riemann 
invariant  will  involve  the  solution  of  an  ordinary  differential  equation  instead  of 
a simple  quadrature  as  was  the  case  in  equation  (A41)  for  Comp  B,  Grade  A.  On 
the  other  hand,  the  case  of  a shock  wave  will  differ  little  in  complexity  from  the 
situation  for  Comp  B,  Grade  A.  The  Sternberg-Walker  equation  is 

p = f1(e)/x  + f2(e)/x3  + f3(e)/x5  + f^(e)/x7  (A47) 

where  the  fj,  i=l,***,4  are  polynomials.  In  this  case,  it  may  even  be  necessary 
to  use  numerical  methods  to  eliminate  e in  order  to  obtain  the  differential 
equation  for  the  isentrope.  In  any  event,  unless  the  choice  of  the  polynomials 
fj  are  extremely  fortuitous,  the  analogue  of  equation  (A46)  will  be  far  more 
complicated.  Although  we  could  set  down  some  more  equations,  the  foregoing  brief 
analysis  should  ind  ate  the  problems  to  be  solved  in  the  case  of  a general 
equation  of  state. 

One  factor  left  to  be  considered  is  that  an  equation  of  state  is  sometimes 
available  only  in  the  form  of  a graph  or  table.  This  may  be  the  case  for 
accurate  computations  involving  real  air  at  extremely  high  temperatures.  Two 
approaches  are  available  in  this  event  for  constructing  solutions  to  the  Riemann 
problem.  First,  the  data  can  be  fitted  by  an  equation  of  state  in  functional  form 
and  the  problem  can  be  solved  as  before.  Second  the  data  can  be  preprocessed  into 
a format  usable  for  table  lookup  in  the  various  steps  of  the  Godunov  iteration. 
Parenthetically,  we  note  that  this  option  is  also  available  for  a very  complicated 
functional  equation  of  state.  The  choice  between  these  two  alternatives  Is 
clearly  a matter  of  computational  efficiency. 

Finally,  we  point  out  that  the  Riemann  problem  is  not  necessarily  well-posed 
for  an  arbitrary  equation  of  state.  This  is  an  area  of  current  interest  in  the 
mathematical  literature.  For  example,  recall  that  in  the  foregoing  we  have 
occasionally  imposed  various  conditions  such  as  unique  solvability  of  nonlinear 
equations  arising  from  the  equation  of  state  (unique  in  the  physical  sense,  not 
mathematically  unique — that  is  the  solution  must  satisfy  the  equation  and  satisfy 
entropy  and  possibly  other  compatibility  conditions).  It  will  be  of  interest  to 
obtain  theorems  showing  the  Riemann  problem  is  well-posed  for  the  equations  of 
state  that  are  used  in  blast  wave  theory.  It  might  be  of  even  greater  interest  to 
obtain  negative  results.  Since  the  Riemann  problem  is  a reasonable  initial 
configuration,  this  would  show  that  the  equation  of  state  is  unphysical. 
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CONCLUSIONS 

We  have  shown  in  this  report  that  Gliiran's  method  is  applicable  to  a variety 
of  problems  involving  strong  shocks,  material  interfaces,  and  blast  waves. 

A priori,  the  method  has  significant  advantages  over  currently  available  finite 
difference  techniques  for  problems  of  this  type.  Whether  or  not  these  advantages 
outweigh  disadvantages  must  await  further  numerical  tests  and  possibly  theoretical 
algorithm  development,  especially  in  multidimensional  situations.  The  last  section 
outlines  the  directions  that  we  would  like  to  see  this  research  take.  The 
numerical  tests  reported  here  using  strongly  nonuniform  waves  are  very  promising 
and  we  expect  similar  results  for  two-dimensional  analogues.  The  appendix  shows 
that  the  extension  of  the  method  to  a general  equation  of  state  is  straightforward, 
but  (potentially)  very  tedious  from  the  point  of  view  of  developing  an  efficient 
computer  code.  In  summary,  there  is  an  excellent  chance  that  Glimm's  method  will, 
in  the  near  future,  be  able  to  yield  qualitatively  superior  results  and  thereby 
j advance  the  field  of  numerical  analysis  of  blast  wave  effects. 
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